Introduction

In this script we are going to map Myeloid subtypes onto the Visium slides.

Libraries

library(Seurat)
library(ggpubr)
library(cowplot)
library(dplyr)
library(ggplot2)
library(stringr)
library(readr)
library(SPOTlight)
library(SPATA2)

Setting parameters

Loading necessary paths and parameters

source(here::here("misc/paths.R"))
source(here::here("utils/bin.R"))

"{myeloid}/{plt_dir}" %>%
  glue::glue() %>%
  here::here() %>%
  dir.create(path = ,
             showWarnings = FALSE,
             recursive = TRUE)

"{myeloid}/{robj_dir}" %>%
  glue::glue() %>%
  here::here() %>%
  dir.create(path = ,
             showWarnings = FALSE,
             recursive = TRUE)

set.seed(123)

Extract sample id and get Donor ID

# sample_id <- params$sample_id
sample_id <- "esvq52_nluss5"
donor_id <- id_sp_df[id_sp_df$gem_id == sample_id, ] %>% dplyr::pull(donor_id)

Load data

We have 8 different datasets that we are going to analyze separately. The spatial data comes from the script 03-clustering/03-clustering_integration.Rmd while the sc data can be found in Ramon’s scRNAseq analysis: /scratch/devel/rmassoni/tonsil_atlas_private/2-DOWNSTREAM_PROCESSING/results/R_objects/processed_seurat_objects/processed_seurat_objects/tonsil_integrated_with_harmony_scrublet_annotated.rds.

sp_obj <- "misc/{robj_dir}/20220215_tonsil_atlas_spatial_seurat_obj.rds" %>%
  glue::glue() %>%
  here::here() %>%
  readRDS(file = .)

# Load SPOTlight data
spotlight_ls <- "{myeloid}/{robj_dir}/spotlight_ls_myeloid.rds" %>%
  glue::glue() %>% 
  here::here() %>%
  readRDS(file = .)

# Load names df data
nm_df <- "{myeloid}/{robj_dir}/myeloid_nm_df.rds" %>%
  glue::glue() %>% 
  here::here() %>%
  readRDS(file = .)

Add colors to cell types

# (nm_df <- data.frame(col_vec))
# nm_df$plt_nm <- rownames(nm_df)
# nm_df$mod_nm <- stringr::str_replace_all(
#   string = rownames(nm_df),
#   pattern = "[[:punct:]]|[[:space:]]",
#   replacement = ".")

Analysis

Preprocess data

decon_mtrx <- spotlight_ls[[2]]
decon_mtrx <- decon_mtrx[, colnames(decon_mtrx) != "res_ss"]

# Set as 0 cell types predicted to be under 3 % of the spot
decon_mtrx[decon_mtrx < 0.03] <- 0

Change column names

new_cn <- data.frame(mod_nm = colnames(decon_mtrx)) %>%
  dplyr::left_join(nm_df, by = "mod_nm") %>%
  dplyr::mutate(plt_nm = dplyr::if_else(is.na(plt_nm), mod_nm, plt_nm)) %>%
  dplyr::distinct() %>%
  dplyr::pull(plt_nm)

colnames(decon_mtrx) <- new_cn

We are going to add the deconvolution to the Seurat object.

sp_obj@meta.data <- cbind(sp_obj@meta.data, decon_mtrx)

Subset sample of interest

sp_sub <- sp_obj[, sp_obj@meta.data$gem_id == sample_id]
sp_sub@images <- sp_sub@images[Seurat::Images(sp_sub) == sample_id]

Look at SPOTlight results

Check Topic profiles

nmf_mod_ls <- spotlight_ls[[1]]
nmf_mod <- nmf_mod_ls[[1]]

h <- NMF::coef(nmf_mod)
rownames(h) <- paste("Topic", 1:nrow(h), sep = "_")
topic_profile_plts <- SPOTlight::dot_plot_profiles_fun(
  h = h,
  train_cell_clust = nmf_mod_ls[[2]])

topic_profile_plts[[2]] +
  ggplot2::theme(axis.text.x = ggplot2::element_text(angle = 90), 
                 axis.text = ggplot2::element_text(size = 12))

Look at all cells profiles

topic_profile_plts[[1]] +
  ggplot2::theme(
    axis.text.y = ggplot2::element_blank(),
    axis.text.x = ggplot2::element_blank(),
    axis.title = ggplot2::element_blank())

Look at cells topic profile

basis_spotlight <- data.frame(NMF::basis(spotlight_ls[[1]][[1]]))

train_labs <- spotlight_ls[[1]][[2]]
colnames(basis_spotlight) <- unique(stringr::str_wrap(train_labs, width = 30))

basis_spotlight[basis_spotlight < 0.0000001] <- 0

DT::datatable(basis_spotlight, filter = "top")

Cell type location

Look at the location of each cell type in each slice separately

# Iterate over cell types
ct <- colnames(decon_mtrx)

# Iterate over images
lapply(names(sp_obj@images), function(nm) {
  print(nm)
  nm_donor <- id_sp_df %>% dplyr::filter(gem_id == nm) %>% dplyr::pull(donor_id)
  # Iterate over cell types
  ct_plt_ls <- lapply(ct, function(i) {
    tmp_plt <- Seurat::SpatialFeaturePlot(
      object = sp_obj,
      features = i,
      alpha = c(0, 1),
      images = nm) +
      ggplot2::scale_fill_gradientn(
        colors = heat.colors(10, rev = TRUE)) +
      ggplot2::scale_alpha(range = c(0, 1)) +
      ggplot2::labs(title = stringr::str_wrap(string = i,
                                     width = 25),
           fill = "") +
      ggplot2::theme(
        plot.title = ggplot2::element_text(
          hjust = 0.5,
          size = 20,
          face = "bold"))
    
    if (sum(sp_sub@meta.data[, i]) == 0) {
      tmp_plt <- suppressMessages(tmp_plt + ggplot2::scale_alpha(range = c(0,0)))
    }
    
    return(tmp_plt)
  })
  
  (plt_arr <- cowplot::plot_grid(
    plotlist = ct_plt_ls,
    axis = "trbl",
    align = "hv",
    nrow = 6,
    ncol = 5))
  
  "{myeloid}/{plt_dir}/cell_type_location_myeloid_{nm_donor}_new.pdf" %>%
    glue::glue() %>%
    here::here() %>%
    cowplot::save_plot(
      filename = .,
      plot = plt_arr,
      base_height = 30,
      base_width = 25)
  })
## [1] "tarwe1_xott6q"
## [1] "c28w2r_7jne4i"
## [1] "esvq52_nluss5"
## [1] "p7hv1g_tjgmyj"
## [1] "gcyl7c_cec61b"
## [1] "zrt7gl_lhyyar"
## [1] "qvwc8t_2vsr67"
## [1] "exvyh1_66caqq"
## [[1]]
## [1] "/scratch/devel/melosua/phd/projects/BCLLatlas/tonsil_atlas/spatial_transcriptomics/myeloid_integration/2020-09-22/plots_2020-09-22/cell_type_location_myeloid_BCLL-2-T_new.pdf"
## 
## [[2]]
## [1] "/scratch/devel/melosua/phd/projects/BCLLatlas/tonsil_atlas/spatial_transcriptomics/myeloid_integration/2020-09-22/plots_2020-09-22/cell_type_location_myeloid_BCLL-8-T_new.pdf"
## 
## [[3]]
## [1] "/scratch/devel/melosua/phd/projects/BCLLatlas/tonsil_atlas/spatial_transcriptomics/myeloid_integration/2020-09-22/plots_2020-09-22/cell_type_location_myeloid_BCLL-10-T_new.pdf"
## 
## [[4]]
## [1] "/scratch/devel/melosua/phd/projects/BCLLatlas/tonsil_atlas/spatial_transcriptomics/myeloid_integration/2020-09-22/plots_2020-09-22/cell_type_location_myeloid_BCLL-12-T_new.pdf"
## 
## [[5]]
## [1] "/scratch/devel/melosua/phd/projects/BCLLatlas/tonsil_atlas/spatial_transcriptomics/myeloid_integration/2020-09-22/plots_2020-09-22/cell_type_location_myeloid_BCLL-13-T_new.pdf"
## 
## [[6]]
## [1] "/scratch/devel/melosua/phd/projects/BCLLatlas/tonsil_atlas/spatial_transcriptomics/myeloid_integration/2020-09-22/plots_2020-09-22/cell_type_location_myeloid_BCLL-14-T_new.pdf"
## 
## [[7]]
## [1] "/scratch/devel/melosua/phd/projects/BCLLatlas/tonsil_atlas/spatial_transcriptomics/myeloid_integration/2020-09-22/plots_2020-09-22/cell_type_location_myeloid_BCLL-9-T_new.pdf"
## 
## [[8]]
## [1] "/scratch/devel/melosua/phd/projects/BCLLatlas/tonsil_atlas/spatial_transcriptomics/myeloid_integration/2020-09-22/plots_2020-09-22/cell_type_location_myeloid_BCLL-11-T_new.pdf"

Boxplots

Prepare data for boxplots

metadata_long <- sp_sub@meta.data %>% 
  # tidyr::pivot_longer(cols = c("annotation"),
  #                     names_to = "stratification_id",
  #                     values_to = "stratification_val") %>%
  tidyr::pivot_longer(cols = all_of(ct), names_to = "ct_key", values_to = "ct_val") %>%
  # dplyr::left_join(col_df, by = c("ct_key" = "ct_name")) %>%
  dplyr::mutate(ct_val = dplyr::if_else(ct_val > 0.001, ct_val, 0))

Box plot of cell type proportion between stratified regions showing the unadjusted ANOVA Pvalue

keep_ct <- metadata_long %>%
  dplyr::group_by(ct_key) %>%
  dplyr::summarise(prop_sum = sum(ct_val)) %>% 
  dplyr::filter(prop_sum > 0) %>%
  dplyr::pull(ct_key)

(bplt <- metadata_long %>%
  dplyr::filter(ct_key %in% keep_ct) %>%
  dplyr::mutate(
    ct_key = stringr::str_wrap(string = ct_key,
                                   width = 30)) %>%
  ggpubr::ggboxplot(data = .,
                    x = "annotation_20220215",
                    y = "ct_val",
                    facet.by = "ct_key",
                    color = "annotation_20220215",
                    fill = "annotation_20220215",
                    add = "jitter",
                    scales = "free",
                    repel = TRUE,
                    outlier.shape = NA,
                    alpha = 0.6,
                    palette = "Set1",
                    ncol = 6) +
  ggplot2::theme(
    strip.text = ggplot2::element_text(size = 18, face = "bold"),
    axis.text.y = ggplot2::element_text(size = 16),
    axis.title.y = ggplot2::element_text(size = 20),
    # axis.text.x = element_text(size = 12, angle = 90, vjust = 0.5, hjust = 0.5),
    axis.text.x = ggplot2::element_blank(),
    legend.text = ggplot2::element_text(size = 18),
    legend.title = ggplot2::element_blank(),
    strip.background = ggplot2::element_blank()) +
  ggplot2::labs(
    y = "Proportion",
    color = "Regions",
    fill = "Regions"))

# bplt <- bplt +
#   ggpubr::stat_compare_means(method = "anova", size = 6) +
#   ggplot2::scale_y_continuous(
#     expand = expansion(mult = c(0, 0.1)),
#     labels = function(x) sprintf("%.2f", x))

"{myeloid}/{plt_dir}/strat_bplot_{donor_id}.pdf" %>%
  glue::glue() %>%
  here::here() %>%
  cowplot::save_plot(
    filename = .,
    plot = bplt,
    base_height = 20,
    base_width = 25)

Cell type correlation matrix

Slide-specific

lapply(Images(sp_obj), function(i) {
    
    se_sub <- subset(sp_obj, subset = gem_id == i)
    # se_sub
    se_sub@images <- se_sub@images[Seurat::Images(se_sub) == i]
    donor_id <- id_sp_df[id_sp_df$gem_id == sample_id, ] %>%
        dplyr::pull(donor_id)
    
    (cor_mtrx_ct <- SCrafty::correlation_heatmap( 
      se = sp_obj,
      feats = ct,
      assay = "Spatial",
      slot = "data") +
       ggplot2::labs(
         title = glue::glue("{donor_id}: Cell-type correlation matrix")))
    
    "{myeloid}/{plt_dir}/cor-mtrx_cell-type_{donor_id}.pdf" %>%
      glue::glue() %>%
      here::here() %>%
      cowplot::save_plot(
        filename = .,
        plot = cor_mtrx_ct,
        base_height = 9,
        base_width = 10)
})
## [[1]]
## [1] "/scratch/devel/melosua/phd/projects/BCLLatlas/tonsil_atlas/spatial_transcriptomics/myeloid_integration/2020-09-22/plots_2020-09-22/cor-mtrx_cell-type_BCLL-10-T.pdf"
## 
## [[2]]
## [1] "/scratch/devel/melosua/phd/projects/BCLLatlas/tonsil_atlas/spatial_transcriptomics/myeloid_integration/2020-09-22/plots_2020-09-22/cor-mtrx_cell-type_BCLL-10-T.pdf"
## 
## [[3]]
## [1] "/scratch/devel/melosua/phd/projects/BCLLatlas/tonsil_atlas/spatial_transcriptomics/myeloid_integration/2020-09-22/plots_2020-09-22/cor-mtrx_cell-type_BCLL-10-T.pdf"
## 
## [[4]]
## [1] "/scratch/devel/melosua/phd/projects/BCLLatlas/tonsil_atlas/spatial_transcriptomics/myeloid_integration/2020-09-22/plots_2020-09-22/cor-mtrx_cell-type_BCLL-10-T.pdf"
## 
## [[5]]
## [1] "/scratch/devel/melosua/phd/projects/BCLLatlas/tonsil_atlas/spatial_transcriptomics/myeloid_integration/2020-09-22/plots_2020-09-22/cor-mtrx_cell-type_BCLL-10-T.pdf"
## 
## [[6]]
## [1] "/scratch/devel/melosua/phd/projects/BCLLatlas/tonsil_atlas/spatial_transcriptomics/myeloid_integration/2020-09-22/plots_2020-09-22/cor-mtrx_cell-type_BCLL-10-T.pdf"
## 
## [[7]]
## [1] "/scratch/devel/melosua/phd/projects/BCLLatlas/tonsil_atlas/spatial_transcriptomics/myeloid_integration/2020-09-22/plots_2020-09-22/cor-mtrx_cell-type_BCLL-10-T.pdf"
## 
## [[8]]
## [1] "/scratch/devel/melosua/phd/projects/BCLLatlas/tonsil_atlas/spatial_transcriptomics/myeloid_integration/2020-09-22/plots_2020-09-22/cor-mtrx_cell-type_BCLL-10-T.pdf"

Integrated

# se_sub <- subset(merged_se, subset = gem_id == "esvq52_nluss5")
# se_sub
# se_sub@images <- se_sub@images[Seurat::Images(se_sub) == "esvq52_nluss5"]

(cor_mtrx_ct <- SCrafty::correlation_heatmap( 
  se = sp_obj,
  feats = ct,
  assay = "Spatial",
  slot = "data") +
   ggplot2::labs(
     title = "Integrated cell-type correlation matrix"))

"{myeloid}/{plt_dir}/cor-mtrx_cell-type_integrated.pdf" %>%
  glue::glue() %>%
  here::here() %>%
  cowplot::save_plot(
    filename = .,
    plot = cor_mtrx_ct,
    base_height = 9,
    base_width = 10)

Session Info

sessionInfo()
## R version 4.0.1 (2020-06-06)
## Platform: x86_64-conda_cos6-linux-gnu (64-bit)
## Running under: Red Hat Enterprise Linux Server release 6.7 (Santiago)
## 
## Matrix products: default
## BLAS/LAPACK: /scratch/groups/singlecell/software/anaconda3/envs/spatial_r/lib/libopenblasp-r0.3.12.so
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C               LC_TIME=es_ES.UTF-8        LC_COLLATE=en_US.UTF-8     LC_MONETARY=es_ES.UTF-8    LC_MESSAGES=en_US.UTF-8    LC_PAPER=es_ES.UTF-8       LC_NAME=C                  LC_ADDRESS=C               LC_TELEPHONE=C             LC_MEASUREMENT=es_ES.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] parallel  stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] tidyr_1.2.0         tibble_3.1.6        Biobase_2.50.0      BiocGenerics_0.36.1 SPATA2_0.1.0        SPOTlight_0.1.7     readr_2.1.2         stringr_1.4.0       dplyr_1.0.8         cowplot_1.1.1       ggpubr_0.4.0        ggplot2_3.3.5       SeuratObject_4.0.4  Seurat_4.1.0       
## 
## loaded via a namespace (and not attached):
##   [1] utf8_1.2.2                  reticulate_1.24             tidyselect_1.1.2            htmlwidgets_1.5.4           grid_4.0.1                  Rtsne_0.15                  munsell_0.5.0               codetools_0.2-18            ragg_1.2.2                  ica_1.0-2                   DT_0.21                     future_1.24.0               miniUI_0.1.1.1              withr_2.5.0                 spatstat.random_2.1-0       colorspace_2.0-3            highr_0.9                   knitr_1.37                  stats4_4.0.1                SingleCellExperiment_1.12.0 ROCR_1.0-11                 ggsignif_0.6.3              tensor_1.5                  listenv_0.8.0               NMF_0.23.0                  MatrixGenerics_1.2.1        labeling_0.4.2              GenomeInfoDbData_1.2.4      polyclip_1.10-0             bit64_4.0.5                 farver_2.1.0                rprojroot_2.0.2             parallelly_1.30.0           vctrs_0.3.8                 generics_0.1.2              xfun_0.30                   R6_2.5.1                    doParallel_1.0.17           GenomeInfoDb_1.26.7         bitops_1.0-7                spatstat.utils_2.3-0        DelayedArray_0.16.3        
##  [43] assertthat_0.2.1            promises_1.2.0.1            scales_1.1.1                vroom_1.5.7                 gtable_0.3.0                globals_0.14.0              goftest_1.2-3               rlang_1.0.2                 systemfonts_1.0.4           splines_4.0.1               rstatix_0.7.0               lazyeval_0.2.2              spatstat.geom_2.3-2         broom_0.7.12                yaml_2.3.5                  reshape2_1.4.4              abind_1.4-5                 crosstalk_1.2.0             backports_1.4.1             httpuv_1.6.5                tools_4.0.1                 gridBase_0.4-7              ellipsis_0.3.2              spatstat.core_2.4-0         jquerylib_0.1.4             RColorBrewer_1.1-2          ggridges_0.5.3              Rcpp_1.0.8                  plyr_1.8.6                  zlibbioc_1.36.0             purrr_0.3.4                 RCurl_1.98-1.6              rpart_4.1.16                deldir_1.0-6                pbapply_1.5-0               S4Vectors_0.28.1            zoo_1.8-9                   SummarizedExperiment_1.20.0 ggrepel_0.9.1               cluster_2.1.2               here_1.0.1                  magrittr_2.0.2             
##  [85] data.table_1.14.2           scattermore_0.8             lmtest_0.9-39               RANN_2.6.1                  fitdistrplus_1.1-6          matrixStats_0.61.0          hms_1.1.1                   patchwork_1.1.1             mime_0.12                   evaluate_0.15               xtable_1.8-4                IRanges_2.24.1              gridExtra_2.3               compiler_4.0.1              KernSmooth_2.23-20          crayon_1.5.0                htmltools_0.5.2             mgcv_1.8-39                 later_1.3.0                 tzdb_0.2.0                  DBI_1.1.2                   corrplot_0.92               MASS_7.3-55                 Matrix_1.4-0                car_3.0-12                  cli_3.2.0                   SCrafty_0.1.0               igraph_1.2.11               GenomicRanges_1.42.0        pkgconfig_2.0.3             registry_0.5-1              plotly_4.10.0               spatstat.sparse_2.1-0       foreach_1.5.2               ggcorrplot_0.1.3            bslib_0.3.1                 rngtools_1.5.2              pkgmaker_0.32.2             XVector_0.30.0              digest_0.6.29               sctransform_0.3.3           RcppAnnoy_0.0.19           
## [127] spatstat.data_2.1-2         rmarkdown_2.12              leiden_0.3.9                uwot_0.1.11                 shiny_1.7.1                 lifecycle_1.0.1             nlme_3.1-155                jsonlite_1.8.0              carData_3.0-5               viridisLite_0.4.0           fansi_1.0.2                 pillar_1.7.0                lattice_0.20-45             fastmap_1.1.0               httr_1.4.2                  survival_3.3-1              glue_1.6.2                  png_0.1-7                   iterators_1.0.14            bit_4.0.4                   stringi_1.7.6               sass_0.4.0                  textshaping_0.3.6           irlba_2.3.5                 future.apply_1.8.1